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Abstract 

The potential of mean force (PMF) between two nano crystals (NCs) represents an effective interac¬ 
tion potential that can be used to study the assembly of NCs to various superstructures. For a given 
temperature, the effective interaction is obtained best from molecular dynamics simulations. Based on a 
density functional approach, this study proposes three methods of predicting the PMF for any tempera¬ 
ture based on a single molecular dynamics simulation for one temperature. The three methods construct 
the PMF by considering the ligands as an ideal gas, as hard-sphere chains, or as Lennard-Jones inter¬ 
action sites. To apply this methodology, the density of the interaction centers must be extracted from 
the simulation data. For the ideal gas model, a straightforward sampling procedure with a hxed lattice 
in space leads to free energies that are too large in order to consistently explain the simulation data for 
different temperatures. Naive sampling does not account for the small momenta added to the NCs when 
coupled to a thermostat. A method is proposed that corrects for the unphysical steps during the simula¬ 
tion. The ideal gas contribution computed for the corrected density is significantly smaller than the one 
obtained from naive sampling and can thus explain the temperature dependence of the PMF correctly. 
For the hard-sphere chain model, where a weighted density is used, the correction of the particle density 
is not essential. However, the PMF calculated based on the corrected density confirms our approach. All 
three models predict PMF curves in very good agreement with simulation results, but they differ in the 
number of input parameters and the computational effort. Based on the modeling results, we predict the 
existence of an additional attractive force at small distances of the NCs — a depletion force. 

PACS 34.20.Gj; 65.80.-g; 82.60.Qr 

Keywords capped nano crystals; thiol ligands; gold core; potential of mean force; perturbation theory; 
depletion attraction 


1 Introduction 

Nano crystals (NCs) — particles the size of a few nanometers, consisting of many thousand atoms, usually 
capped with ligands — are building blocks of new materials with applications in various fields of modern 
engineering [13]. Much about them, especially about their interaction with other NCs, is known from 
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Figure 1: Two exemplar NCs and their PMF. Panel A illustrates the considered system of two Aui 47 (SCg )58 
NCs at distance R. The NCs have a gold core consisting of 147 Au atoms (forming an icosahedra, yellow), 
58 alkanethiol ligands consisting of one SH- (brown) and eight CH^j-groups (cyan). Panel B shows the PMFs 
for five temperatures as obtained from MD simulations as well as the statistical errors. 


experiments [28] and molecular simulations [20], including simplified models [15]. Due to their complex 
structure and interaction behavior, a good theoretical knowledge is often missing. The precise shape of pair- 
and many-body-potentials, which determine the assembly of NCs into superstructures, is not well understood 
theoretically and is so far obtained best by molecular simulations [12,21]. 

In this article we investigate the effective pair potential employing models of density functional theory 
(DFT) while relying on data from molecular dynamics (MD) simulations. The procedure allows for an 
accurate prediction of effective pair potentials at temperatures other than the temperature considered for 
the MD simulation. To assess our approach, we run simulations with two thiol capped gold NCs for various 
temperatures (Fig. 1). Due to the large number of interaction sites (1-lOk) and the requirement of sufficiently 
long simulation runs computations of this kind are expensive. Therefore, we try to reduce the number of 
necessary simulations and use all the information that is generated — not only the final estimates of the 
mean force. 

For instance, we analyze trajectory data from our MD simulations, determine the density of the relevant 
interaction sites at the ligand caps and apply their density to various fluid models of density functional 
theory (DFT). In doing so we neglect some of the NCs’ structure. For example, considering a chain model, 
we neglect the adsorption of the thiol ligand to the gold core, and considering an ideal gas, we even neglect 
molecular forces. How can such an approach be justified? 

The potential of mean force (PMF) will be represented by the difference of two Helmholtz free energies 
[11]; the minuend is the free energy of the two NCs at focal distance, the subtrahend is the free energy of the 
NCs at infinite distance. Therefore, any interaction that does not vary in strength when changing the NCs’ 
distance to each other will not contribute to the PMF. One example are the forces by which the ligands are 
absorbed to the gold surface, another example are the chemical bonds between the atoms in the thiols. That 
is, when calculating the PMF, one can neglect most of the NCs’ structure. 

2 Methods 

2.1 MD simulations 

For a set of five temperatures T, constrained MD simulations have been performed to determine the PMF 
^t{R) between two exemplar NCs (Fig. 1). Following a united-atom approach with the Lennard-Jones force 
helds from Schapotschnikow et al. [20], we consider the SH-, CH 2 - and CHg-groups in the thiol ligands as sin¬ 
gle interaction sites. Simulations were conducted as described in our previous work [5], using the GROMACS 


2 

















software [23] with a leap-frog stochastic dynamics integrator of step size At = 2 ps and friction constant 
7 = .5ps“^. This integrator also acts as a thermostat. For every temperature (T = 300K,..., 500K), 10 
runs were performed in the NVT ensemble, each with an equilibration time of 1 ns and a total runtime of 
t = 10 ns, providing us with trajectories of 5000 coordinate snapshots. 

The centers of mass of the two NCs were aligned with the a;-axis of the simulation box and their distance 
ri 2 was fixed during a simulation by freezing the gold atoms in space. Only in between simulations the 
distance was changed, incrementally (Ari 2 = .5nm), from large = 5nm > dc + 2L) to small distances 
= 2 nm > dc)\ c?c = 1.8 nm is the diameter of the core and L = 1.1 nm is the length of the ligands. Here 
we make two implicit assumptions [5]: (i) the rotation of the NCs is negligible [20] and (ii) the core-to-core 
interaction is unimportant compared to the ligand-to-ligand and ligand-to-core interactions [25]. 

To obtain the PMF, we evaluated the forces, Fi and F 2 , that act on the centers of mass of the gold cores 
over the course of a simulation run. The mean force between the two NCs is given by 

= , (1) 

where the angular bracket denotes the average in the constrained canonical ensemble and Bx the unit vector 
in x-direction. The potential at distance ri 2 = R results from integration over larger distances, 

noo 

$T(i?)= / Fm{r)dr, (2) 

Jr 

which in our discrete setting translates into ^t{R) = X]ra=o Ari 2 , with Fm{'r'^ 2 ^) = 

0. Consequently, statistical errors become larger towards smaller distances R (cf. Fig. IB). 

2.2 Thermodynamics 

As a special case of Jarzynski’s non-equilibrium equality [11], the PMF obtained from our MD simulations 
can be represented by the difference of two Helmholtz free energies, 

$T(-R) = A(i?)-A(oo), (3) 

describing thermodynamic states of the NCs separated by a certain distance R [5]. In the context of the 
simulations, the limit A(oo) in Eq. (3) corresponds to a value at finite distance, A{r^^). 

Most naturally [32], the Helmholtz energy is expanded in powers of the inverse temperature /3 = \/{kBT), 

00 

^A = ^/3"a„, (4) 

n—0 

where ks defines Boltzmann’s constant. Even if not indicated, the coefficients a„ also depend on /3, but in 
the considered temperature range the dependence is unimportant. 

The coefficients a„ can be calculated using theoretical methods of statistical physics (see Sect. 2.5 below). 
In practice, very few of them need to be known to reproduce the experimental data. For dispersive molecular 
interactions (e.g., the Lennard-Jones interactions between the molecular groups of our NCs), the first two 
are sufficient (Fig. 2AB), 

/3A = ao + Pai+Oif3^). (5) 

Without experimental or theoretical knowledge, simulations for two [5] or more temperatures /?„ = 
l/^ksTn) are necessary to estimate the two coefficients. By applying the pseudo-inverse p = {mA ■ m)~^ ■ mA 
of the Vandermonde matrix m = ((1, /3i), (1, /32), • ■ ■) , 

( 00 , 01 )’’’ = p ■ {I3 iAt-^,P2At2,- , (6) 


3 






Figure 2: Prediction of the expansion coefficients. Panels A and B show the expansion coefficients (dots) 
obtained from MD simulation data of five temperatures via the pseudo-inverse. These coefhcients fit the 
PMFs of the MD simulation best (C). The area (red) around the best fit contains the coefficients obtained 
for pairs (gray) of temperatures; the border lines represent the envelope of all ten pairs of temperatures 
and thus indicate the uncertainty for our PMF calculation. Panel D shows the fit for just two PMF curves 
(300 K, 500 K); this is the result of [5]. The “prediction error” represents the mean square deviation (first 
value); the second value represents the error restricted to NC-distances i? > = 2.35nm. 


we obtain a fit of the coefhcients for any given number of simulations while ensuring that the Htting error 
e = (3A — oo — j3ax is minimal, • £:| —min (cf. Fig. 2C). 

For di-polar interactions, second order expansion is needed, etc., and the number of necessary simulations 
and coefhcients a„ in Eq. Eq. (5) increases. Again, the best ht is achieved through Eq. Eq. (6) in a higher 
dimension. 


2.3 Interpolation and extrapolation 

Going back to the linear expansion Eq. (5), let us consider the case of two unknown coefhcients. Then the 
pseudo-inverse reduces to the usual inverse. When inserting the resulting coefhcients (cf. (2.14-15) in [5]) 
into the expansion. 


/3$j’ = 


(/3 — I32)Pi ^Ti + (/3i — /3)/32‘hT2 

Pi - P2 


(7) 


one identihes the solution as linear interpolation for /3 $t. Furthermore, one realizes that in accordance 
with the high temperature expansion the interpolation parameter is given by the inverse temperature. The 
interpolation results are shown in Figure 2D. 
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Figure 3: PCA-correction. Panel A shows the three eigenvalues of the inertia tensor (for a representative 
NC-distance of i? = 3nm) after each time step of the MD simulation. The values fluctuate around one 
small and two (very similar) large numbers. Panel B shows the eigenvectors that correspond the largest 
eigenvalues over the simulation time; the eigenvectors to the smallest eigenvalue have been aligned to the 
x'-axis (the common x-direction); cylindrical symmetry is recognizable. Panel C illustrates the alignment 
of the lattice in accordance with the center of mass and the smallest eigenvector of the inertia tensor. This 
procedure must be applied to the trajectory data after each time step of the MD simulation. 


Now let us assume that one coefficient is given. Then only one simulation (performed at temperature Tq, 
say) is needed to calculate If o-o is known, 


= {13 — Po) Aoo + $To ) 


( 8 ) 


and if oi is known, 

Aoi + ^ $To ■ (9) 

These formulas are readily derived from Eq. (5); note that the coefficients require the following adjustment, 
Aci^(./^) — ci^{R{ o.,2(c5o). 

2.4 Physical particle density 

To reduce the number of simulations, we will calculate the expansion coefficients in Eq. (5) using three 
models of density functional theory (Sect. 3). By applying this methodology we must determine the density 
of the relevant interaction sites. Naive sampling does not account for the momenta that are added to the 
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nano system when coupled to a thermostat. During each time step At, the stochastic dynamics integrator 
adds randomly chosen momenta to the interaction sites [26] . 

The center of mass and the tensor of inertia calculated from our simulation data fluctuate notably 
(cf. Fig. 3A). Therefore we propose a procedure that corrects for the unphysical interference (cf. Fig. 3C). 
We re-adjust the center of mass and the tensor of inertia of our trajectory data by applying principal 
component analysis (PCA) after each time step of a simulation run. That is, we move the lattice (Fig. 3C) 
so that effectively we sample a density (referred to as physical) that keeps the center of mass and the tensor 
of inertia (i.e., because of rotational symmetry only the lowest eigenvalue) of the nano system fixed. Here 
averages over the 5000 time steps of a simulation run serve as reference values. In our simulations the gold 
cores are kept frozen in space, therefore we only correct the data of the ligand caps. 

Performing PCA (using the LAPACK software package [1]), we observe one small- and two very similar 
large eigenvalues (Fig. 3A). The eigenvectors of the two large eigenvalues form an almost steadily distributed 
circle around the x-axis (cf. Fig. 3B). Therefore it is reasonable to propose cylindrical symmetry for the 
density of the interaction sites in the cap, even though the cores of the NCs are icosahedras [29]. This 
provides another, practical argument why performing our simulation by freezing the NCs’ cores in space 
(and not only restraining their centers of mass and allowing for rotations) is an acceptable simplification [5]. 

2.5 Perturbation theory 

The coefficients of the high temperature expansion Eq. (4) are determined by the molecular interaction. We 
recall that the partition function Q = N\)~^ J exp {—defines the Helmholtz free energy 

A = -p-HnQ ( 10 ) 

in terms of the molecular potentials ^ denotes the number of interaction sites, 

z)(®'l^(-) the pair potential between sites i and j, and A the thermodynamic wavelength. Here, perturbation 
theory, which model the transition of the free energy from a well-known reference fluid (with pair potential 
Vo) to a target fluid (with pair potential v), defines the coefficients a„, n > 0 [10]. The high-temperature 
expansion, Eq. (4), reads 


(3A = ao + l3{W)+0{l3^), (11) 

where W = represents the perturbation and {W) = f W exp(—/3V[)^°^) dr^/ J exp(—/3V[}°^) dr^ 

denotes the corresponding average for the canonical ensemble of the reference fluid. Examples are content 
of the subsequent section. 

3 Results 

Three thermodynamic models will be used to calculate one of the expansion coefficients in Eq. (4). Relying 
on linear expansion Eq. (5), MD simulations for only one reference temperature (e.g., Tg = 350 K) are needed 
to determine the PME for any other temperature T; see Eq. (8) and Eq. (9). We reconstruct PMFs for the 
available set of simulations and explain the intricacies when applying the models. We start off with the 
simplest one. 

3.1 The ideal gas model 

Here we assume that the entropic contribution (i.e., the coefficient ag = —S/ks) represents an ideal gas. It 
is, however, the ideal gas contribution of a heterogeneous fluid with a density distribution obtained from the 
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MD simulation. Given a continuous particle density, p = lim , as /S.V 0, the ideal gas contribution is 
known (for classical fluids) as 

ao = / p{r) (In {A^p{r)) - l) cfr . (12) 

This formula can be simplified when only taking differences of the Helmholtz free energy, 

ao = y p{r)liip{r)(fr . (13) 

According to Eq. (3), the difference concerns two states: one, where the NCs are located at a distance R 
apart from each other and, another, where one of two NCs is separated to infinity. We thereby assume that 
the intra-molecular partition sum is unchanged for the two states modeled with Eq. (13). 

The result represents Shannon’s entropy [18]. Its only input is the particle density, which can be obtained 
from trajectory data. While considering differences of the free energy, we restrict ourselves to interaction 
sites that markedly change the position with respect to (the center of mass of) their NC when moving the 
NCs in space. These interaction sites are the CHa;-groups (i.e., CH 2 - and CHa-groups) of the thiol chains. 

The simplicity of this approach is deceiving. The density in our formulas is a continuous quantity, but the 
positional coordinates of the interaction sites obtained from an MD simulation are determined for discrete 
time steps. We find that a cubic lattice with mesh size S = .05 nm allows for sampling a reliable discrete 
density, leading to reproducible values for the Helmholtz free energy. Yet these values turn out to be ca. 20% 
too large compared to the predicted values (cf. Fig. 4A). 

We believe the discrepancy is due to the thermostat. The stochastic dynamics integrator used in the 
simulations randomly adds momenta to the molecular groups. This is exactly what we observe from the 
trajectory data. The center of mass and the tensor of inertia of the interaction sites fluctuate over time 
(Fig. 3AB), despite the fact that the molecular groups of the two NCs are the only particles simulated in 
space. 

To correct for the unaccounted physical interaction we utilize a density obtained through PCA, ending 
up with an entropic contribution that is significantly smaller than the one obtained from naive sampling. 
Above a certain distance of the NCs, i? > = 2.35 nm, the agreement of our prediction with the oo-data 

(as determined from MD simulations at various temperatures) is very good (Fig. 4A). 

In order to appreciate this result, we note that an accurate prediction of Uq (as shown in Fig. 4A) leads 
to a reliable estimate of the PMF for any temperature (via Eq. (8)) using MD simulation data of only one 
temperature. The predicted PMF at various temperatures are shown in Fig. 4B. 


3.2 Reference fluid of hard-sphere chains 


A much more advanced model for the thiol ligands of our NCs is given by Wertheim’s theory of hard-sphere 
chains [31]. Here the entropic contribution represents a chain with m hard-sphere segments, containing 
Carnahan-Starling’s term [6]. In contrast to the previous approach (ideal gas model), the ideal gas term is 
neglected [8] here, because now the reference ideal gas state is one where interaction sites are considered free 
(unconstrained to chain structure). The ideal gas contribution then cancels in the difference of Eq. (3). The 
reference fluid (term of zero order) is 


ao 



/ 477 — 377^ 


(1 — 1 / 777 ) In 


1-^/2 \ 

(1-77)3; 


d^r 


(14) 


where t] = VaP denotes the packing fraction and a the hard-sphere diameter. According to Tarazona [24], a 
weighted particle density 


p{r) 


= J p{r')0{cF 


\r — r'j) dr' 


(15) 
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Figure 4: Comparison of PMFs from MD simulations and PMFs predicted from three models: the ideal gas 
(AB), the hard-sphere chains (CD), and the Lennard-Jones fluid (EF). Panels AC and E show the coefficients 
ao{R) and ai{R), resp. Panels BDF display the PMF curves predicted from the models, that are represented 
by blue lines in the corresponding diagrams ACE. For the ideal gas method (A), the coefficients relying on 
a naive particle density (green, dotted) are up to 20% larger than the coefficients based on the physical 
particle density (blue). Above a critical distance (i?° = 2.35nm), the blue coefficient curves approximate 
the expected coefficients (red dots) very well. Below the critical distance the deviation becomes large. For 
the two other methods, the critical distance is higher {Rl = 2.7nm). For the Lennard-Jones fluid, BH’s 
approach (blue) performs better than WCA’s (green, dotted). For all methods, the reference temperature 
To = 350 K is used. Fluctuations in the calculated coefficient curves (< 5%) were smoothed out. 

















is utilized, computed from the physical density p{r). 

To determine the two parameters, m and cr, we use a homo-segmented group contribution method [19], 
where the thiol chains are composed by a fractional number (m = 3.804, by Eq. 2 in [19]) of associating hard- 
sphere segments with equal diameter (cr = .386nm, by Eq. 3 in [19]) and hence equal volume {V^ = ttct^/G). 
Figure 4CD shows the results. Similar to the ideal gas model, the fit is very good for distances of the two 
NCs above = 2.7 nm. 

Lo Verso et al. [15] also used a chain model for a similar but artificial nano system, where the density p{r) 
was sampled naively. In contrast to our ideal gas model, no involvement of the thermostat and no dependence 
on the mesh size has been reported. In fact, this coincides with our observations. When applying the naive 
density to our chain model we obtain almost the same results as when applying the physical density. This 
could be due to the weighted particle density and its smoothing behavior. One has to be careful though 
when implementing Eq. (15) on a lattice; the lattice must be fine enough to sufficiently represent a single 
hard-sphere volume Vo- by summation over its cells. 

3.3 Lennard-Jones fluid 

In our third model we take a different approach. Rather than estimating the reference fluid oq, we now 
estimate the ai term (via Eq. (9)). The approach is similar in practical application, because the respective 
other term (oq or oi) is determined by MD simulations and Eq. (5). We consider the ligand segments 
as a Lennard-Jones fluid. That is, we apply the force field that was used in the MD simulation, v{r) = 
4e ((cr/r)^^ — {a/r)^). However, we only consider the CHa;-groups as interaction sites, and we only compute 
the interaction term in first order perturbation theory, 

Oi = i JJ p{r)p{r')g{\r-r'\)vi{\r-r'\)d^r(fr', (16) 

where vi = v — vq denotes the perturbation and g{-) the pair distribution of the reference fluid. As reference 
fluid we consider hard-spheres, and numerically we approximate their pair distribution using the Percus- 
Yevik’s method [22]. We apply two standard theories: the one by Barker-Henderson (BH) [3] and the one 
by Week-Chandler-Andersen (WCA) [30]. The density p{r) is sampled again from the simulated trajectories 
after PCA is performed. The results are shown in Figure 4EF. The fitting is very good above the distance 
Rl — 2.7 nm. Naive sampling produces slightly worse results. 

BH performs better than WCA. This is a bit surprising, as BH’s higher order terms — unlike WCA 
— usually also contribute to Eq. (4). Maybe the choice of the reference fluid with an only temperature 
dependent hard-sphere diameter is responsible. BH’s reference potential is steeper than WCA’s and thus 
closer to a hard-sphere fluid. 

The computational effort of summing over all interactions Eq. (16) is much higher than in the previous 
two methods. We took advantage of cylindrical and mirror symmetry when sampling the particle density, 
reducing the computation time by the factor 1/64. 

4 Discussion 

For an exemplar system of two thiol capped gold core NCs we estimated the PMF using three models. 
Above a critical distance (i? > i?/), each of the three models led to PMF curves of satisfactory quality. The 
methods underlying these models differ considerably with respect to the number of input parameters and 
the computational effort. 


9 



B 


rs 


rc 


0.5 


1.0 1.5 

r / nm 


2.0 



Figure 5: Radial density. Density profiles of NC 1 und NC2 practically coincide. In Panel A, the NCs are 
separated by a distance of i? = 5mn. For smaller separation distances (cf. B) the densities have shifted 
— away from intermediate radial distances, for the SH-group, and towards larger distances, for the first 
CH 2 -group; that is, SH makes room for secondary CH 2 -groups of the other NC. 

4.1 Model comparison 

The ideal gas model returned the best PMFs and requires no extra parameters. Yet one must correct the 
simulated particle distribution via PCA to obtain a physical density. As MD simulations usually involve 
(virtual) interaction with a thermostat, our correction should be relevant for many similar problems. One 
disadvantage of this methodology is the additional computational cost (when performing PCA for 5000 data 
sets), but it is far less than the cost of possibly additional MD simulations. 

The hard-sphere chain model requires knowledge of two additional parameters, the hard-sphere diameter 
and the chain length. Both these parameters were obtained independently from a homo-segmented group 
contribution method. The fit of the PMF-data is almost as good as with the ideal gas model. For further 
improvement, one may try to employ a hetero-segmented method [19]. Because of the weighted density 
that enters the formulas, the model turned out to be robust with respect to small imprecisions of the 
particle density. That is, compared to the ideal gas model, no PCA correction is needed. Therefore, the 
computational cost is smaller too. 

The ’Lennard-Jones model’, did not perform as well as the other two models, in particular for small 
NC-distances. There are several possible reasons. The radial pair distribution function for spheres is only 
a rough approximation for the pair distribution of heavily structured NCs [25]. For densely overlapping 
ligands (i.e., NC-distances R < 2.7nm), the sum over the Lennard-Jones interactions turned out too high. 
Furthermore, we introduced a lattice of rather big mesh size, because summation over all the attractive 
forces is costly and extends with cubic power. The proposed rotational and mirror symmetries for the lattice 
entries will have contributed to errors as well. 

In each of the three models, the NC ligand layer is studied as a fluid. Taking this aspect of the approach 
seriously and regarding the ligand layer as a solvent for the core, we are able to clarify the role of an additional 
force observed at small distances of the NCs, as explained now. 


4.2 Model extension: depletion attraction 

The change in monotony of the entropic coefficient ao(i?) around the critical distance R = and the 
relatively small error bounds (cf. Fig. lA) suggest that there is an additional attractive force Fq{R) not yet 
accounted for in our models. Let us assume that this force behaves according to Hook’s law with spring 
constant k — with onset at distance Rq > i?[! and becoming larger at smaller distances of the two NCs. The 


10 
























R/nm R/nm 


Figure 6: Effect of depletion forces on the PMF. Expansion coefficients and PMFs of the ideal gas model and 
the Lennard-Jones model are shown with and without contribution due to depletion. The onset distance is 
estimated by theoretical arguments and thus chosen to be the same for both model fluids, Rq = 2.357 nm. 
This value is very close to the best fit obtained for the ideal gas. The parameter k is fitted to each model, 
individually. 
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additional contribution to the PMF would then be quadratic, 

=-'^{Ro - Rf 0{Ro - R), (17) 

with 9 being the Heaviside function. We adjusted the parameters, Rq and k, to the PMF-data generated by 
the MD simulation. The best result is achieved for a onset distance Rq = 2.374 nm slightly above the critical 
distance = 2.35 nm; the regression shows excellent agreement with the data (Fig. 6A), with error values 
similar to [5]. The harmonic ansatz function, Eq. (17), represents the simplest form of a depletion force [2,9]. 
In practice, adjusting the parameter kappa requires one additional MD simulation. The additional simulation 
is in conflict with our objective to reduce the simulation effort to a single MD run, and for many applications 
it is sufficient to simply neglect the short-ranged potential, Eq. (17). Nonetheless, from a more fundamental 
viewpoint we are interested in elaborating the force further, showing that the short-ranged force contribution 
is indeed likely the result of a depletion force. In doing so, we show that the distance parameter Rq can be 
estimated from geometric data. 

Depletion forces between two NCs are apparently attractive forces. The effect is caused by the depletion 
of smaller entities (here interaction sites of the ligands) between the NCs due to geometrical size exclusion. 
As a result, the local density of ligands between the NCs is taking on small values, leading to a low local 
pressure. The depletion forces are thus attractive. Depletion has been used to explain phenomena in the 
context of nano-sized particles — with polymer parts taking over the role of solvent particles [7, 14,16,17]. 

In our nano system, only the CH^j-groups move and thus mimic solvent particles. Therefore, depletion 
attraction is supposed start at (or slightly above) the distance where only one CH 2 -group fits in between the 
two NCs. 

The considered CH 2 -group is bound to a SH-group, which is bonded rather immobile to the gold core 
of one NC (#1, say). The surface of the other gold core (of NC2) is densely covered with SH-groups 
(cf. Fig. 8B). The average distance of the two NCs, i?sc)(Sj where the first CH 2 -group of NC 1 touches the 
SH-group surface of NC2, should then be an estimate for the onset distance Rq. In fact, we found that 
numerically the optimal onset distance coincides with 


77sc)(S —’’C + <5c)(s + — 2.357nm (18) 

up to 17pm (which is far below mesh size 6 = 50pm); here rs = .950nm and rc = 1.107nm are the mean 
radii of the SH- and the first CH 2 -groups about the gold cores (obtained from trajectory data; cf. Fig. 5), 
and <5c)(s = -312 nm represents the average distance between the centers of CH 2 (NC 1) and SH (NC2) over 
the bumpy SH-surface of NC2 (calculated in Appendix A). 

4.3 Outlook 

To test the range of validity we performed MD-simulations for pairs of NCs with other core sizes and ligand 
lengths, and we also studied heterogeneous NC pairs [5]. For all these cases the modeling approach is in 
good agreement with a series of independent MD results (cf. Fig. 7). It needs to be examined, however, if for 
more extreme dimensions (e.g., very long ligands) and for other core and/or cap materials our methodology 
leads to comparably good results. 

After two-body interaction, three-body interaction must be studied [4] . 


A Depletion: the onset distance 

Looking at the geometric structure of the NCs, one cannot recognize a simple rule that determines the 
positions of the SH-groups adsorbed to the Au-icosahedra. It is known, however, that for the NCs considered 
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Figure 7: Modeling results for another pair of NCs, Aui 4 i 5 (SCi 2 ) 242 - The NCs consist of gold cores with 1415 
Au atoms and 242 alkanethiol ligands (one SH- and twelve CHa;-groups). Panel A shows the coefficient ao(i?) 
based on MD data for three temperatures; the best estimate is indicated by red dots, the corresponding error 
margin by the surrounding red area. Panel B compares the predicted PMF (lines) and the PMF obtained 
from MD simulations (dots). Panels C and D show the results of the ideal gas model with (and without) 
depletion; Tq = 400 K. The onset distance [Rq = 4.284 nm) has been chosen according to our geometric 
analysis, based on the mean radii (rs = 1.894 nm, rc = 2.052 nm) obtained from the trajectory data. Here, 
the onset distance causing the best ht (4.445 nm) is notably higher and likely to be attributed to the longer 
ligands. 
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Figure 8: Contact distance between two NCs. The blue surface z{x, y) in Panel A illustrates the center of the 
first CH 2 -group of NC 1 when touching (at least) one of three SH-groups of NC2 with centers at zero-level 
(red triangle). Panel B shows the geometric parameters involved. 

here, there are H = 58 ligands adsorbed to the gold surface and that this number represents the densest 
packing [20]. For the analysis in this appendix, we consider Lennard-Jones interaction sites as hard spheres 
of diameter a. 

For simplicity, we assume that the SH-surface is a large sphere on which the SH-groups are regularly 
distributed. That is, one would propose that the SH-surface is tiled by hexagons with a SH-group in the 
center. Then one only needs to study three SH-groups, and the only unknown parameter is the distance e 
between the boundaries of two neighboring SH-spheres (cf. Fig. 8). 

Assuming that the hexagonal tiling covers all the SH-groups distributed around the gold core at mean 
distance rs, 

4TTrl=HA, (19) 

where A = ( °’sh-i-£ ^^ represents the hexagon area surrounding one SH-group, the distance between the 

boundaries can be estimated by 

Regarding depletion, we propose that the first CH 2 -groups (following SH) are the parts of the ligands 
that undergo depletion for small NC-distances. Over a triangular region, as illustrated in Figure 8, the 
hard-sphere diameters ((Tch 2 = -396nm and ctsh = .445nm [20]) determine the contact distances between 
the first CH 2 -group of NC 1 (cyan) and three neighboring SH-groups of NC 2 (red). These contact distances 
(i.e., before depletion starts) are calculated by 

z{x, y) = + _ ^)2 ^ (21) 

This formula only applies to a particular part of the equilateral triangle. Due to symmetry, it is sufficient 
to average over one sixth of the surface (e.g., the triangle given by the points (0, Ay), (0, 0), (Ax, 0) with 
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Ay = (crs + £^)/2 and Ax = Ay/i/S), 


2 pAx pip(x) 

<5c)(S = y J z{x,y)dydx - Sc - Ss = .300 nm (22) 

(where •(/'(x) = Ay — ^ x defines the hypotenuse of that triangle; illustrated by dots in Fig. 8A). Note that 
the center of the three SH-spheres is slightly below their mean radial distance from the gold core; the offset 
is given by (5s = rg — \Jr^ — Ax'^ — Ay^ = .054 nm. A similar correction is included for the CH 2 -sphere, 
Sc = re - r{x,y)dydx = .011 nm, where locally r{x,y) = - {Ax - x)'^ - y^ defines 

the projection of rc onto the core-to-core direction. 
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